<N 

o 

Oh 
<D 



0) 

i 

a 

a 
o 
o 



> 

r-- 
m 

m 

en 

o 



X 



Quasiparticle spectra from a non-empirical optimally-tuned 
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We present a method for obtaining outer valence quasiparticle excitation energies from a DFT-based calcula- 
tion, with accuracy that is comparable to that of many-body perturbation theory within the GW approximation. 
The approach uses a range- separated hybrid density functional, with asymptotically exact and short-range frac- 
tional Fock exchange. The functional contains two parameters - the range separation and the short-range Fock 
fraction. Both are determined non- empiric ally, per system, based on satisfaction of exact physical constraints for 
the ionization potential and many-electron self-interaction, respectively. The accuracy of the method is demon- 
strated on four important benchmark organic molecules: perylene, pentacene, 3,4,9, 10-perylene-tetracarboxylic- 
dianydride (PTCDA) and 1,4,5,8-naphthalene-tetracarboxylic dianhydride (NTCDA). We envision that for finite 
systems the approach could provide an inexpensive alternative to GW, opening the door to the study of presently 
out of reach large-scale systems. 



Development of a non-empirical theory for quantitative 
electronic structure calculations, which combines predictive 
power with computational simplicity, is a long-standing chal- 
lenge for molecular and solid-state physics dEJ. Presently, 
many -body perturbation theory within the GW approximation 
BH3 is widely considered to be the first principles approach 
that provides the best balance between accuracy and computa- 
tional tractability. This approach is couched in a formally rig- 
orous theory for quasiparticle excitations and has been shown 
to provide remarkably quantitative predictions for the elec- 
tronic structure of a wide variety of molecular, solid-state, and 
low-dimensional systems (see, e.g., BHTOl). 

Unfortunately, present-day GW calculations are still signif- 
icantly limited in system size and complexity. They can also 
be challenging to converge l lTl - fl3ll . Therefore, it is common 
practice to rely instead on density functional theory (DFT) 
|[T4ll . which is much simpler computationally. However, this 
comes at a significant cost in accuracy. Solutions of the Kohn- 
Sham equation (in either its original [15] or generalized lfT6l 
form) generally do not rigorously correspond to quasiparticle 
energies and orbitals. Practical DFT calculations can still be, 
and often are, successful because occupied DFT eigenvalues 
can, in principle, serve as good approximations to removal 
energies of energetically high-lying occupied orbitals (51 Wl\ - 
l20l . Even so, two major problems remain ifTTl . First, it is 
often found that the eigenvalue spectrum can depend strongly, 
and even qualitatively, on the choice of the approximate den- 
sity functional. Second, the energies of the highest occupied 
molecular orbital (HOMO) and the lowest unoccupied molec- 
ular orbital (LUMO) typically do not correspond to the ion- 
ization potential and electron affinity, respectively. 

In this Letter, we show that DFT-based calculations, in 
which outer- valence orbitals do represent quasiparticle exci- 



tations, are in fact possible, opening the door to inexpensive 
prediction of quasiparticle excitations. Our approach is based 
on a range-separated hybrid density functional, which is op- 
timally tuned to obey Koopmans' (ionization potential) the- 
orem and to minimize many-electron self-interaction errors, 
without any recourse to empiricism. 

Recently, Stein et al. (21] suggested a new method for pre- 
dicting the fundamental gap of finite systems from general- 
ized Kohn-Sham HOMO and LUMO eigenvalues, based on 
a non-empirical optimally-tuned range- separated hybrid (OT- 
RSH) functional. In an RSH functional, the Coulomb repul- 
sion is partitioned into a short-range (SR) and a long-range 
(LR) part, such that the LR exchange is treated with a Fock 
operator whereas the SR exchange is treated using (semi- 
)local exchange (22). The range- separation parameter, 7, is 
optimally-tuned by demanding that the DFT version of Koop- 
mans' theorem (ionization potential theorem) be obeyed, i.e., 
by determining 7, per-system, such that the HOMO eigenval- 
ues of the neutral and anionic system are as close as possible 
to the ionization potential and electron affinity of the neutral, 
respectively, making the gap contribution of a derivative dis- 
continuity (T71 [23l [24) negligible. Refaely-Abramson et al. 
have shown the efficacy of this approach, using the optimally- 
tuned RSH, denoted here as OT-7RSH, for a range of organic 
molecules of relevance to photovoltaics (25ll . This results in 
HOMO and LUMO levels on par with the GW ones. It can 
perhaps be hoped, then, that the outer-valence occupied or- 
bitals of finite systems could also be well described with this 
tuned functional. 

To test this notion, Fig. [TJa,b) presents a comparison of 
experimental photoemission spectra and theoretical DFT and 
GW eigenvalue spectra [ 26 ] for perylene and pentacene. Here 
and throughout, DFT eigenvalue spectra were performed us- 



ing QChem [27] with the cc-PVTZ basis set [28) and GW 
calculations were performed using the Berkeley GW code 
[29] (details are given in the supplementary information, SI). 
For both molecules, the GW spectrum agrees well with ex- 
periment. The DFT spectrum constructed using the gener- 
alized gradient approximation (GGA) in its Perdew-Burke- 
Ernzerhof (PBE) form [ 30 ] does not agree with GW, but is pri- 
marily rigidly shifted from it, owing to the missing derivative- 
discontinuity in the exchange-correlation functional [23]. Hy- 
brid functionals, e.g., PBE0[31 ] (which is based on the PBE 
functional but with 25% of the GGA exchange replaced by 
Fock exchange) mitigate (but do not solve) the derivative dis- 
continuity problem [17, 24] and consequently improve the 
spectrum by shifting (and slightly stretching) it. Still, the 
PBEO spectrum exhibits a significant rigid shift. Conversely, 
the OT-7RSH is in quantitative agreement with GW and ex- 
periment (average unsigned error of ~0.2 eV over a range of 
-3 eV below the HOMO), as hoped for. 

Unfortunately, this simple idea is not sufficient for 
more complex molecules, such as 3,4,9, 10-perylene- 
tetracarboxylic-dianydride (PTCDA) and 1,4,5,8- 
naphthalene-tetracarboxylic dianhydride (NTCDA) 

(Fig. [TJc,d)), where GW spectra agree with experiment 
but the OT-7RSH spectra present serious deviations from 
GW in both orbital position and ordering. These molecules 
have been chosen because, while still reasonably simple, 
they exhibit a mixture of localized (on the anhydride side 
groups) and delocalized (on the perylene or naphthalene 
core) outer valence orbitals, as shown in Fig. [TJc,d). Dori 
et al. [32] pointed out that this causes the PBE spectrum of 
PTCDA to be in poor agreement with GW even after a rigid 
shift, as also shown in Fig. [TJc). The spectral distortions 
result primarily from spurious positive energy shifts of the 
localized orbitals, leading to the conjecture that they reflect 
significant self-interaction errors (SIE)[32|. Korzdorfer et al. 
l33l showed that NTCDA has a similar problem, as can be 
seen in Fig. |TJd), and proved the conjecture by quantifying 
the per-orbital SIE for both molecules. Furthermore, they 
showed that self-interaction-corrected calculations, within a 
generalized optimized effective potential scheme, provide a 
non-empirical route for obtaining agreement with experiment, 
up to a rigid shift of both HOMO and LUMO (and possibly 
some mild stretching). 

For both PTCDA E2 and NTCDA (33, a different non- 
empirical route for improvement of the eigenvalue spectral 
shape is the use of the above-mentioned PBEO hybrid func- 
tional, which possesses a fixed fraction of exact exchange, as 
shown in Fig. [IJc,d). Indeed, Korzdorfer and Kiimmel l34l 
showed that such hybrids can incorporate an important part of 
the first order correction of the Kohn-Sham eigenvalues when 
used in a generalized Kohn-Sham way, i.e., with non-local 
Fock exchange. However, recent orbital tomography experi- 
ments [ 35 ] showed that the orbital ordering in hybrids can still 
be wrong. This is also reflected in Fig.[T]:, and is likely a con- 
sequence of the hybrid functional not being self-interaction 
free. In addition, conventional hybrid functionals do not re- 



solve the ionization potential and fundamental gap problem. 
Thus, our goal is to provide a generalized Kohn-Sham scheme 
that does yield the correct excitation thresholds and is also 
self-interaction free. We now show that this is an achievable 
goal. 

The results of Fig.[T]suggest that further improvement could 
be obtained by combining a fraction of SR Fock exchange that 
would improve the description of the localized orbitals, with 
LR Fock exchange that is essential for gap prediction. Such 
generalization of the RSH scheme has in fact been suggested 
by Yanai et al. (37), who partitioned the Coulomb operator 
according to 



1 a + /3erf (77*) 
r r 



1 - [a + /3erf (77*)] 



(1) 



Note that this form reduces to that of a conventional hybrid 
functional with the choice f3 = (for PBEO, a = 0.25) and to 
a "pure LR" RSH for a = 0, f3 = 1. 

The parameters a, (3, and 7 can be determined as universal 
parameters semi-empirically, as done for all three of them in 
the CAM-B3LYP functional of Yanai et al. E3 However, as 
discussed in, e.g., refs. EQ El E3 [38ti40l. no set of fixed 
values is universally useful for spectroscopy. Instead, we pur- 
sue the optimal tuning strategy, where we determine all three 
parameters from first principles, per-system, based on satis- 
faction of physical constraints. 

First, following ref. [41], we insist on a + f3 = 1. This 
choice guarantees that full Fock exchange is obtained asymp- 
totically, which enforces the correct asymptotic potential. 
This, in turn, is essential for retaining accurate gap predictions 
(24). As in Refs. (42l|43), we shall use PBE-based semi-local 
exchange and correlation components. However, in these arti- 
cles a was taken as a constant of and 0.2, respectively, with 
7 a universal empirically determined constant. We shall seek 
to optimize both a and 7 non-empirically, based on additional 
constraints. 

Per each choice of a, 7 can be determined from first princi- 
ples by enforcing the ionization potential (Koopmans') theo- 
rem, i.e., by choosing 7 such that the HOMO eigenvalue is as 
close as possible to the ionization potential [44, 45]. In princi- 
ple, this exact condition should be obeyed for any stable ionic 
state of the molecule. Therefore, one can seek 7 that best 
satisfies (say, in the least squares sense) multiple ionization 
potential conditions, by minimizing a target function J(j) of 
the form: 



J 2 (7; a) 



Z^( £ fi(N+i) 



■IP^ a (N + i)) 2 



(2) 



where en(N-\-i) is the HOMO eigenvalue of the N+i electrons 
system, N being the number of electrons in the neutral sys- 
tem and i an integer representing electrons added or removed 
from it, with IP(N + i) the corresponding ionization poten- 
tial, calculated from energy differences. In previous work that 
emphasized accurate gap prediction |2T][24l|25l[38j|40l, 7 was 
chosen so as to satisfy this condition as closely as possible for 
both neutral and anion (i.e., i =0 and 1 in Eq. [2]), so as to 
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FIG. 1. (Color online) Outer valence eigenvalue spectra of (a) Perylene, (b) Pentacene, (c) PTCDA, (d) NTCDA, as obtained from DFT 
calculations using the PBE, PBEO, and OT-7RSH, with additional OT-o^RSH results for PTCDA and NTCDA, compared with GW eigenvalue 
spectra and with experimental gas phase photoemission spectra |36|. H and L denote HOMO and LUMO, respectively. For PTCDA and 
NTCDA, representative localized and delocalized orbitals are presented. All computational spectra have been broadened by convolution with 
a gaussian to facilitate comparison with experiment. 



obtain both the ionization potential and the electron affinity 
(the latter being equal to the ionization potential of the anion). 
Here, this does not suffice, as 7 must also reflect a balance of 
SR and LR exchange appropriate for the treatment of localized 
states. Because with OT-7RSH the highest localized orbital is 
HOMO-1, for both PTCDA and NTCDA, we additionally im- 
pose an ionization potential condition for the cation, i.e., z=-l, 
0, and 1 in Eq. [2] (see SI for additional details). 

The remaining question, then, is how to determine the SR 
Fock exchange fraction, a. To understand the effect of a on 
the spectrum, Fig. [2] shows the outer valence eigenvalues as a 
function of a, for the example of PTCDA (similar results for 
all other molecules are given in the SI). For each choice of a, 
the optimal value of 7, determined by employing Eq. [2] with 
a triple summation, has been used and is also shown. Sev- 
eral important trends can be distinguished immediately. First, 
as a increases, the optimized 7 decreases. This is reason- 
able: the range above which the exchange is dominated by its 
LR contribution roughly corresponds to I/7, and the extent 
of LR Fock corrections should decrease with increasing SR 
Fock contributions. Second, for a between and 0.5, 7 tuning 
is successful throughout in maintaining a HOMO-LUMO gap 
that is constant to within ~0.05 eV and is in excellent agree- 
ment with GW. Larger a values are not given because for too 
large a determining a corresponding 7 that obeys Koopmans' 
theorem to a meaningful accuracy is no longer possible. This 
makes sense, because there is a limit to the extent of SR Fock 
exchange that can be used while still maintaining compatibil- 
ity with a semi-local correlation expression tV7\ . 

Third, Fig. [2] clearly exposes the different behavior of 
the two types of orbitals present in the outer valence re- 
gion. Eigenvalues corresponding to delocalized orbitals (on 



the perylene core) are essentially indifferent to the choice of a 
(to within a mean value of 0.2 eV). Conversely, all anhydride- 
localized orbitals are highly sensitive to a. As an example, 
the doubly-degenerate orbital, that is HOMO- 1/2 for a=0, 
is HOMO-5/6 for a=0.5, and it changes in energy from ~- 
1.5 eV to ~-2.8 eV. A similar picture emerges for NTCDA. 
For perylene and pentacene, all outer valence orbitals within 
~3 eV below the HOMO are delocalized and the spectrum is 
largely independent of a (see SI for details), which explains 
the success of the SR-exchange-free OT-7RSH functional for 
these molecules (Fig. l(a,b)). Deeper lying orbitals possess 
different degrees of localization and are outside the scope of 
this work. 

How to choose an optimal a, then, without empiricism? 
Recently, Srebro and Autschbach l46l suggested that it can 
be obtained by insisting on an additional property satisfied by 
exact DFT, namely, that for ensemble states described by a 
fractional number of electrons the total energy versus particle 
number curve must be piecewise linear [44]. They used such 
tuning to obtain an accurate CuCl electric field gradient. We 
stipulate that satisfaction of the piecewise linearity constraint 
is important for spectroscopy as well, for two reasons. First, 
Yang and co-workers have emphasized the importance of lin- 
ear segments for accurate gap prediction (47] |48j . Second, 
enforcing piecewise linearity has been shown to be essential 
for the accurate spectroscopy of localized states |49l [50) and 
deviation from this condition is in fact often dubbed a "many- 
electron SIE" or a "derealization error" l48ll5Tl . 

Several groups have already shown that for a well- 
constructed RSH functional, curves of the energy as a func- 
tional of the fractional number of electrons, for the [N-1,N] 
and [N, N+l] segments, are much more linear than those ob- 
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FIG. 2. (Color online) Eigenvalue energy as a function of the short- 
range Fock fraction, a, with the optimal value of the range- separation 
parameter, 7 (Bohr -1 ) deduced for each choice of a, for PTCDA. 
HOMO and LUMO levels are marked. Rectangle denotes the opti- 
mal a, determined from the minimization procedure described in the 
text. 



tained with conventional functionals, even in the absence of 
optimal tuning l24l[39ll46ll5Tti53l . We have performed frac- 
tional electron calculations for our benchmark molecules us- 
ing NWChem [ 54 ] with the same basis set as above. As shown 
in Fig. fusing PTCDA as an example, the above findings ap- 
ply here as well: Whereas PBE and PBEO exhibit a notable 
deviation from linearity, for an OT-RSH functional the devi- 
ation from linearity is too small to be detected by the naked 
eye. An alternative approach to assessing segment linearity, 
which is more directly relevant to spectroscopy, is to consider 
the dependence of the eigenvalues on the fractional number 
of electrons iBTlL as shown in the inset of Fig. [3] for [N-2, 
N+l]. If the linear- segment constraint is satisfied, the HOMO 
eigenvalue should be constant between integer electron val- 
ues, owing to Janak's theorem l55l . Again, only the OT-RSH 
functional obeys this requirement closely enough. The opti- 
mal value of a can thus be obtained by choosing the a (and 
therefore the corresponding 7) that minimizes the three cur- 
vatures of the AE(AN) curve for -2 < A7V < 1. It is found 
to be 0.2 for both PTCDA and NTCDA. 

Satisfactorily, we find that for both molecules the orbital or- 
dering for the optimal a is similar to the GW ordering (with 
the exception of PTCDA HOMO-5 orbital which is slightly 
misplaced), without any need for level shifting: OT-a^RSH 
eigenvalues of PTCDA delocalized states deviate from the 
same in GW by ~0.18 eV (the largest deviation being 0.25 
eV). For localized states, this deviation is ~0.01 eV. For 
NTCDA, the average deviation of all states from GW is ~0.07 
eV. These numbers are well within the accepted accuracy of 
either calculation. Indeed, the OT-a^RSH spectrum, also 
shown in Fig. [T] agrees extremely well with the GW one for 
both absolute HOMO and LUMO positions and the quasipar- 
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FIG. 3. (Color online) Deviation of total energy from that of the 
neutral molecule, AE, and HOMO eigenvalue, shomo (inset), as a 
function of the fractional deviation of the number of electrons from 
that of the neutral molecule, AN, computed for PTCDA using PBE, 
PBEO, and OT-cryRSH (a=0.2, 7=0.160 Bohr" 1 ). The table shows 
the curvature of each functional, in eV, obtained from fitting the 
AE(AN) curve to a second order polynomial. 



tide spectrum of filled states, for all examined molecules. 

In conclusion, we demonstrated, using four important 
benchmark molecules, that DFT-based calculations can reach 
an accuracy that is comparable to that of GW calculations. 
This was achieved by using a PBE-based range- separated hy- 
brid density functional, with asymptotically exact and short- 
range fractional Fock exchange. Importantly, both range- 
separation and Fock fraction are determined non-empirically, 
based on satisfaction of exact constraints for the ionization po- 
tential and many-electron self-interaction error, respectively, 
resulting in full predictive power for the outer valence elec- 
tronic structure. We envision that the approach could be use- 
ful directly as a low-cost alternative to GW that offers good 
accuracy for both frontier and non-frontier quasiparticle ex- 
citation energies. Additionally, because perturbative "one- 
shot" Go Wo is known to be sensitive to the DFT starting point 
O [56] [53, our approach provides a novel optimal starting 
point for subsequent GW calculations. 
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